##

rm(list = ls())

## 

library(tidyverse)

## Load saved results

res_df_rdd_muni <- read_rds('data/rdd_municipal_results.rds') %>% 
  mutate(unit = 'Municipal election') %>% 
  filter(str_detect(outcome, 'agg|turnout'))

res_df_rdd_bundestag <- read_rds('data/rdd_federal_results.rds') %>% 
  mutate(unit = 'Federal election') %>% 
  filter(str_detect(outcome, 'agg|turnout')) %>% 
  mutate(outcome = str_remove(outcome, "_party"))

res_df_rdd_landtag <- read_rds('data/rdd_state_results.rds') %>% 
  mutate(unit = 'State election') %>% 
  filter(str_detect(outcome, 'agg|turnout'))

## Combine and do some cleanup / renaming

results_rdd <- list(res_df_rdd_landtag, 
                    res_df_rdd_bundestag, 
                    res_df_rdd_muni)%>% 
  reduce(bind_rows) %>%
  filter(!str_detect(outcome, 'turnout')) %>% 
  mutate(analysis = 'RDD') %>% 
  mutate(bw_type = ifelse(bw %% 500 == 0, 'regular', 'opt_bw')) %>%
  filter(!str_detect(outcome, 'cand')) %>% 
  mutate(exclude_state = ifelse(exclude_state == "State excluded",
                                'B-W excluded', 'Full sample')) %>% 
  mutate(outcome= recode(outcome, 
                         `agg_left` = 'Left-wing\nparties',
                         `agg_center` = 'Centrist\nparties',
                         `agg_right` = 'Right-wing\nparties')) %>% 
  mutate(unit = factor(unit, levels = unique(unit)[c(2,1,3)])) %>% 
  mutate(outcome = factor(outcome, levels = unique(outcome)[c(1,3,2)]))

# Figure A.6: Main results ----

pd <- position_dodge(0.5)

p3 <- ggplot(results_rdd %>% 
               filter(str_detect(exclude_state, 'B-W')),
             aes(outcome, estimate)) +
  geom_hline(yintercept = 0, linetype = 'dotted',
             color = 'grey40')  +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0, size = 0.6) +
  geom_errorbar(aes(ymin = conf.low90, ymax = conf.high90), 
                width = 0, size = 1) +
  geom_point(shape =21,
             fill = 'white', size = 2.5) +
  xlab('') + ylab('RD estimate\n(percentage points)') + theme_bw() +
  facet_wrap(~unit, ncol = 3) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))+
  scale_y_continuous(breaks = seq(-10, 5, 2.5))
p3

## Get results for different bandwidths 

res_bw <- read_rds('data/rdd_federal_results_sensitivity.rds') %>% 
  filter(!str_detect(outcome, 'turnout')) %>% 
  mutate(bw_type = ifelse(bw %% 500 == 0, 'regular', 'opt_bw')) %>%
  mutate(exclude_state = ifelse(exclude_state == "State excluded",
                                'B-W excluded', 'Full sample')) %>% 
  mutate(outcome= recode(outcome, 
                         `agg_left_party` = 'Left-wing parties',
                         `agg_center_party` = 'Centrist parties',
                         `agg_right_party` = 'Right-wing parties')) %>% 
  mutate(outcome = factor(outcome, levels = unique(outcome)[c(1,3,2)])) %>% 
  filter(bw < 5001) %>% 
  filter(str_detect(outcome, "Left|Centrist|Right"))

## Subset

diff_bw <- results_rdd %>% 
  filter(str_detect(unit, 'Federal')) %>% 
  mutate(outcome= recode(outcome, 
                         `Left-wing\nparties` = 'Left-wing parties',
                         `Centrist\nparties` = 'Centrist parties',
                         `Right-wing\nparties` = 'Right-wing parties'))


## Combine w/ main results using optimal BW

plot_df <- bind_rows(res_bw, diff_bw) %>% 
  mutate(outcome = factor(outcome, levels = unique(outcome)[c(1,2,3)]))

# Figure A.9 : BW sensitivity

p1 <- ggplot(plot_df, aes(bw, estimate, bw_type)) +
  geom_hline(yintercept = 0, linetype = 'dotted',
             color = 'grey40')  +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high,
                    color = bw_type), position = pd,
                width = 0, size = 0.65) +
  geom_errorbar(aes(ymin = conf.low90, ymax = conf.high90,
                    color = bw_type), position = pd,
                width = 0, size = 1.05) +
  geom_point(position = pd, size = 2.5,
             shape = 21,
             aes(fill = bw_type)) +
  scale_color_manual(values = c('grey50', 'black')[2:1],
                     labels = c('Alternative bandwidth',
                                'MSE-optimal bandwidth')[2:1],
                     name = '') +
  scale_fill_manual(values = c('white', 'black')[2:1],
                    labels = c('Alternative bandwidth',
                               'MSE-optimal bandwidth')[2:1],
                    name = '') +
  xlab('') + ylab('RD estimate (percentage points)') + theme_bw() +
  facet_grid(exclude_state~outcome) +
  scale_x_continuous(breaks = seq(0, 10000, 1000))+
  theme(legend.position = 'bottom') +
  xlab('Bandwidth') + ylab('RD estimate (percentage points)') +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
p1